#ifndef FIELD_INCLUDED
#define FIELD_INCLUDED 1

#include "types.h"

typedef struct
{
	int        nx;
	int        ny;
	double   **array;
	const char *name;
} Field;

#define CHECK_FIELD(f) \
assert((f)!=NULL); \
assert((f)->nx>0); \
assert((f)->ny>0); \
assert((f)->array!=NULL); \
assert((f)->array[0]!=NULL); \
assert((f)->name != NULL)

#define CHECK_FIELDS(f,g)     \
assert( (f) != NULL );        \
assert( (g) != NULL );        \
assert( (f)->nx == (g)->nx ); \
assert( (f)->ny == (g)->ny )

#define ARRAY_OF(f) ((f)->array)
#define CONST_ARRAY_OF(f) (const double **) ARRAY_OF(f)

//! create a field
Field *field_create(const char *name,
                    const int   nx,
                    const int   ny );

//! create a field with the same size
Field *field_align( const char *name, const Field *f);

//! delete a field
void field_delete(Field *f);

//! set all values of the field
void field_set(Field        *f,
               const double  value);

//! display field info
void field_info(const Field *f,
                FILE        *fp);

//! save to vtk
void field_write(const Field  *f,
                 const Field  *g,
				 const Field  *h,
				 const Field  *sig11,
				 const Field  *sig12,
				 const Field  *sig22,
				 const Field  *P,
                 const Field  *dvg,
				 const Field  *eta1, 
				 const Field  *zeta1,
				 const int cycle,
                 const double  dx,
                 const double  dy);

//! save to vti, version 2
/*void field_save(const Field *f,
                const char  *prefix,
                const int    cycle,
                const double dx,
                const double dy);*/


//! compute laplacian with PBC
void field_laplacian(Field        *f_laplacian,
                     const   Field *f_source,
                     const double   inv_dx2,
                     const double   inv_dy2
                     );

//! compute laplacian with rho_u/rho
void field_laplacian_rhu(Field      *f_laplacian_rhu, 
						 const  Field *f_source,
						 const  Field *f_rho,
						 const double inv_dx2, 
						 const double inv_dy2
						 );
//! compute laplacian with rho_v/rho
void field_laplacian_rhv(Field     *f_laplacian_rhv,
						 const Field *f_source, 
						 const Field *f_rho, 
						 const double inv_dx2, 
						 const double inv_dy2
						 ); 
//! compute first part of non-linear equation with rho_u_u and rho_u_v
void field_nonlinear_1(Field   *f_rx, 
					   const Field *f_source1,
					   const Field *f_source2, 
					   const Field *f_rho, 
					   const double inv_dx, 
					   const double inv_dy
					   );
//! compute second part of non-linear equation with rho_v_v and rho_u_v
void field_nonlinear_2(Field  *f_ry, 
					   const Field *f_source1, 
					   const Field *f_source2, 
					   const Field *f_rho, 
					   const double inv_dx, 
					   const double inv_dy
					   );  

//! calcul de la tension de surface et aussi la partie de la viscosite.
void field_tension(Field *f_sigma_1_1,
				   Field *f_sigma_1_2,
				   Field *f_sigma_2_2, 
				   const Field * f_rho1,
				   const Field * f_rhu1,
				   const Field * f_rhv1,
				   const Field * f_eta1,
				   const Field * f_zeta1,
				   const double K, 
				   const double inv_dx, 
				   const double inv_dy,
				   const double inv_dx2, 
				   const double inv_dy2); 


//! calcul de la tension de surface 
void field_tension_bis(Field *f_sigma_1_1,
					 Field *f_sigma_1_2,
					 Field *f_sigma_2_2, 
					 const Field * f_rho1,
					 const double K, 
					 const double inv_dx, 
					 const double inv_dy,
					 const double inv_dx2, 
					 const double inv_dy2); 


//! compute gradient with PBC
void field_gradient(Field *      f_grad_x,
                    Field *      f_grad_y,
                    const Field *f,
                    const double inv_dx,
                    const double inv_dy);

//! compute divergence with PBC
void field_divergence(Field       *f_dvg,
                      const Field *f_x,
                      const Field *f_y,
                      const double inv_dx,
                      const double inv_dy);
//! calcul de la divergence d'un tenseur. 
void field_divergence_tenseur1(Field *f_dvg1,
							   const Field *f_x,
							   const Field *f_y,
							   const double inv_dx,
							   const double inv_dy);
void field_divergence_tenseur2(Field *f_dvg2,
							   const Field *f_x,
							   const Field *f_y,
							   const double inv_dx,
							   const double inv_dy);

//! set pressure according to the state equation
void field_set_pressure(Field          *f_P,
                        state_equation  eq,
                        const Field    *f_rho,
                        const double    T);

//!< on calcule la viscosite en fonction de la densite
void field_set_viscosity(Field         *f_eta1, 
						 state_equation visco, 
						 const Field   *f_rho, 
						 const double T);

//!< on calcule la deuxieme viscosite en fonction de la densite
void field_set_second_viscosity(Field *f_zeta1,
								state_equation second_visco,
								const Field *f_rho, 
								const double T);

//! compute momenta from rho and velocities
void field_compute_momenta(Field       *f_Ju,
                           Field       *f_Jv,
                           const Field *f_rho,
                           const Field *f_u,
                           const Field *f_v);

//! compute velocities from rho and momenta
void field_compute_velocities(Field       *f_u,
                              Field       *f_v,
                              const Field *f_rho,
                              const Field *f_Ju,
                              const Field *f_Jv);
								 
//! on calcule l'increment.
void field_compute_increment(Field *f_delta_Ju,
							 const Field *f_dvg,
							 const Field *f_grad,
							 const Field *f_rx);

void field_compute_increment_1(Field *f_delta_Ju1,
							   const Field *f_dvg, 
							   const Field *f_grad,
							   const Field *f_laplacian,
							   const Field *f_rx,
							   const double eta); 

//! add field contents
void field_add( Field *f_a, const Field *f_b );

//! sub field contents
void field_sub( Field *f_a, const Field *f_b, const double lam );

//! mul/add field contents
void field_muladd( Field *f_a, const double lam, const double eta, const Field *f_b, const Field * f_c, const Field *f_d, const Field *f_dvg);
void field_muladd_Ju(Field *f_a, const double lam, const Field *f_b);


#endif

